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Abstract. The dispersionless Kadomtsev-Petviashvili (dKP) equation {ut + 
uux)x = Uyy is one of the simplest nonlinear wave equations describing two- 
dimensional shocks. To solve the dKP equation we use a coordinate transformation 
inspired by the method of characteristics for the one-dimensional Hopf equation 
Ut + uux = 0. We show numerically that the solutions to the transformed equation 
develops singularities at later times with respect to the solution of the dKP equation. 
This permits us to extend the dKP solution as the graph of a multivalued function 
beyond the critical time when the gradients blow up. This overturned solution is 
multivalued in a lip shape region in the (x, y) plane, where the solution of the dKP 
equation exists in a weak sense only, and a shock front develops. A local expansion 
reveals the universal scaling structure of the shock, which after a suitable change 
of coordinates corresponds to a generic cusp catastrophe. We provide a heuristic 
derivation of the shock front position near the critical point for the solution of the 
dKP equation, and study the solution of the dKP equation when a small amount 
of dissipation is added. Using multiple-scale analysis, we show that in the limit of 
small dissipation and near the critical point of the dKP solution, the solution of the 
dissipative dKP equation converges to a Pearcey integral. We test and illustrate our 
results by detailed comparisons with numerical simulations of both the regularized 
equation, the dKP equation, and the asymptotic description given in terms of the 
Pearcey integral. 
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1. Introduction 


Perhaps the best known example of a singularity in an evolution equation is the 
formation of jump discontinuities of the density and of the velocity field in the Euler 
equations of compressible gas dynamics. As these discontinuities propagate, they are 
known as shock waves. In the case of a planar shock front, the problem can be reduced 
to a one-dimensional equation for the velocity alone [HO] (a so-called simple wave). The 
resulting wave profile overturns to form an s-shaped curve, the point where the gradient 
first becomes infinite (known as the gradient catastrophe) corresponds to the formation 
of a shock. From the overturned solution the physical solution can be reconstructed 
by inserting a jump discontinuity (the shock). The shock solution is a weak solution of 
the equation, which satisfies additional conditions motivated by physical considerations 
EH. This shock solution is also found by taking the limit of vanishing viscosity in the 
dissipative form of the equations, yielding a weak solution (see jU] for conservation 
laws in one space dimension and m [19] for hyperbolic equations in several space 
dimensions). 

The existence of such gradient catastrophe points has been proved in [21 ES] for 
hyperbolic equations in many space dimensions. However, to the best of our knowledge, 
if the initial condition depends on two or three spatial variables, little is known about 
the two or three-dimensional spatial structure of the shock near the blow-up points of 
the gradients. In particular, it would be interesting to know the self-similar structure 
of the solution both before and after shock formation [Tl| . A rare instance of where we 
have a more or less complete understanding of a higher dimensional singularity is the 
spatial structure of caustics of wave fronts in the approximation of geometrical optics 
mm- Two-dimensional wave breaking has also been studied in |13|, using a simple 
kinematic equation, for which an exact implicit solution is available. 

In this paper, we study the formation of two-dimensional shocks in a simple 
nonlinear wave equation known variously as the dispersionless Kadomtsev-Petviashvili 
(dKP) equation | 22 |, or the Zabolotskaya-Khokhlov (ZK) equation [121 . The equation 
has the advantage that its one-dimensional form, the Hopf equation, has only one 
family of characteristics. The dKP equation can be seen as a long wavelength version 
of the original Kadomtsev-Petviashvili (KP) equation [22] : 

{ut + uu^ -I- "^XXX^X iUyjy, (1-1) 

but with the highest order dispersive term u^xx dropped, namely 


(u^ T ^^yy- 

The subscript denotes the derivative with respect to the variable. With a -|- sign on 
the right hand side, (1.1) is known as the KPI equation, or as the KPII equation 
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in the opposite case. However, in the case of the dKP equation the two signs are 
equivalent under the transformation u —>■ —u and x —>■ —x, and for the remainder of 
this manuscript we will consider only the positive sign. Depending on context, the 
KP equation describes wave profiles for layers of inviscid fluid of finite depth, waves in 
plasmas, or the propagation of sound beams in nonlinear media. 

While the Cauchy problem for the KP equation is globally well-posed in a suitable 
space [lO], by dropping the dispersive term, the dKP equation becomes a nonlocal 
scalar conservation law in two space variables. Even for smooth initial data, the 
solution remains smooth only for finite time. In [IH] it is shown that the solution 
of the dKP equation is locally well posed in the Sobolev space H^, s > 2, so that for 
s > 4 one has classical solutions. Particular solutions of the dKP equations have been 
obtained with several techniques naiiziETiiisiiis]. The Cauchy problem for the dKP 
equation and shock formation have been studied recently in [31ISSIISHIEZI, using the 
inverse scattering transformation, which relies on the integrability of the dKP inherited 
from the KP equation gZllsol. 

To sketch a derivation of the dKP equation, we follow the original derivation of 
the KP equation [22] . We start from the Hopf equation 




UUx = 0 


( 1 . 2 ) 


for a wave field u, with only a convective non-linearity. This is the simplest model 
equation describing wave steepening and shock formation. In a frame of reference 
moving at the sound speed c, a simple wave can be shown to be described by (1.2) [3(1] . 


Assuming a weak y-dependence, we add a small correction ip on the right hand side of 


( 1 . 2 ); 


Ut + UU^ = Ip. 


(1.3) 


For a wave of small amplitude, the second term in the above equation can be neglected. 
Assuming a dispersion relation uj = kc = y/k^ + kyC, one obtains in a frame of reference 


moving along the x-axis with velocity c that u = kc — k^c ~ cky/{2kx). For (1.3) to 
match this dispersion relation, we must have ip^ ~ cUyy/2. Taking the ^-derivative on 


both sides of (1.3) we obtain 


T uUx)x 2^^yy 

Rescaling x —>■ —x u —>■ —u and y —)■ one arrives at the equation 


{lit T 'W^xfx 


yy 


(1.4) 


Note that in spite of its name, the dKP equation (1.4) contains dispersion, and only 


the highest order dispersive term has been dropped relative to fill. Other contexts 


in which (1.4) is used are described in |H]. 
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The Hopf equation ( |1.2[ ) is solved by observing that the velocity is constant along 
characteristic curves given by [9]: 




(1.5) 


Thus for any initial condition u{x, 0) = Uo{x), one hnds an exact solution u{x, t) = Uo(^) 
in implicit form. Wave breaking occurs when two characteristics cross, which always 
occurs when the initial condition has negative slope. A shock hrst forms along the 
characteristic originating from the point of greatest negative slope by absolute value, 
where the solution u(x, t) has a point of blow up of the gradient. 

Thus if one expands the initial condition about one finds that the profile 
assumes a characteristic s-shape m 

Ax — AuAt + t\u'''{^f)Avf’/Q = 0, (1.6) 


where Au = u — Uc, and Ax = x — Xc — u^it — tf). For At = t — tc > 0 (after shock 
formation), the profile has become multivalued. Balancing the three terms in ( |1.6[ ), 
one sees directly that Au must be of order and so Ax of oder [131 El. 


If one solves (1.4) with an initial condition which depends on y, the equation can 


no longer be solved with the method of characteristics. The idea underlying this paper 
is that the dependence on the ^/-coordinate is weak, so the structure of the solution is 
essentially the same as before, but the different stages of overturning are “unfolded” 
in the y-direction [33] • This means that effectively the singularity time becomes a 
function of y. If we choose the origin such that a singularity occurs at y = 0 first, and 
expand tc in a Taylor series near y = 0, we obtain tc{y) — tc(0) + -|- 0{y^), with 

a > 0 a constant and 4 ( 0 ) = tc- This means that At = t — tc — ay"^ = f— ay"^, and the 
two-dimensional wave breaking is governed by the scalings 


Au ~ Ax ~ Ay ~ 


(1.7) 


In this paper we will show that the scalings (1.7) indeed describe the similarity structure 
of wave breaking in the dKP equation. 


The estimates (1.7) imply that Ay ^ Ax near the shock, consistent with our 


assumption of a slow variation in the ^-direction. The central idea of our paper is to 


use this insight to generalize the characteristic transformation (1.5) to allow for a slow 
y-dependence: 

u{x,y,t) = F(^,y,t) 

X = tF{^,y,t)+^ 


( 1 . 8 ) 


Applying transformation (|1.8) to (1.4[) results in a PDF for F(^,y,t) which we will 


study in the next section (see equation (2.8)); the initial condition for F is given by 


F{x,y,0) = UQ{x,y). 


(1.9) 
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Note that if the initial data Uo{x,y) has no y-dependence, (1.8) yields the exact 
characteristic solution with F{x,y,t) = uo{x) as described before; in particular, F 
is y and time-independent. As in the method of characteristics, the solution u{x,y,t) 
of the dKP equation encounters a gradient catastrophe when the transformation x = 
tF{^, y,t) + ^ defining ^ = i{x, y, t) is not invertible, namely when tF^{^, y,t) + 1 = 0. 
Our numerical results show that as a result of the unfolding (1.8), the function F{^,y,t) 
remains regular at the time tc of shock formation of the solution u{x, y, t) of the dKP 
equation. Moreover, our numerics indicate the derivatives of F remain bounded for 
times substantially beyond tc- However, since F satisfies a nonlinear equation (see (2.8) 
below), we believe that F will typically develop a singularity for some time t > tc, we 
give an example of such a singularity in a particular case. 


Manakov and Santini [35|, [38] have proposed a transformation for analysing the 
gradient catastrophe of dKP equation which is superficially similar to ours, which is 
motivated by the inverse scattering transform. Their transformation differs from ours 
by a factor of 2 in front of the unfolding term: 


{ u{x,y,t) = F{C,y,t) 

x^=2tFic,y,t) + C (1.10) 

F{C,y,0) = uo{C,y) 


as a result, the transformation does not unfold the overturned profile if there is no 
^-dependence. In fact, transformation of the Hopf equation leads to the same equation 
Ft — FFi^ = 0 as before, but with propagation in the opposite direction, and with 
the same initial data F{f,0) = uo(C). This means that for y-independent initial data 
localized in the x-direction, F(f, t) will experience a gradient catastrophe before u(x, t) 
does, if the initial profile is steeper on the left than on the right. The same remains 
true for solutions of the full dKP equation with localized initial data: we checked 
numerically that for the initial data considered in this manuscript, i.e. the x-derivative 
of a Schwartz function, the function F{(, y, t) suffers a gradient catastrophe before a 
gradient catastrophe occurs in the original profile u(x,y,t). 


To further illustrate the difference between the two parameterizations, note that 
combining (1.8) and (1.10) one finds F in terms of F: 


I F(ly,t)=F{(,y,t) 

\ f = tF(c.y,t)+C 


or F in terms of F\ 


F{C,y,t) = F{^,y,t) 
C = -tF(^,y,t) +^. 


( 1 . 12 ) 


If we assume that F{(,y,t) has no singularities and that 2tF(^{f,y,t) -|- 1 > 0 in some 
time interval [0,t'], then it follows from (1.10) that the solution u{x,y,t) of the dKP 
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equation is regular in the same time interval. But since we also have tFi^{f, y,t) + l > 0, 
it follows from (1.11) that is regular in [0, t'] as well. 

On the other hand, assuming that F{f^,y,t) is regular and tF^{f^,y,t) + 1 > 0 in 
some time interval [0,t'], it follows from (1.8) that once again u{x,y,t) is regular in 
[0, t']. However, this does not imply that F{f, y, t) is regular, since it may happen that 
—tF^(^, y,t) + l = 0 for some t G (0, t'], even thoug h tF^ {^, y,t) + l > 0 for all t G [0, t']. 
This argument shows that F{(,y,t), as defined by (1.10), might encounter singularities 
even before u{x, y, t) does. 


Our formulation allows us to find spectrally accurate solutions to F = F(^,y,t), 
from which u{x,y,t) can easily be reconstructed. The alternative would be to 
use numerical methods for hyperbolic equations which remain stable even after the 
formation of shocks [ 32 ] • However, these methods introduce numerical dissipation near 
the shock, which renders the solution inaccurate. These sources of inaccuracy can be 
avoided using our transformation. The main results of this paper are the following: 


• in section 2 we describe the solution of the dKP equation by using a transformation 
inspired by the method of characteristics and by |35| . This transformation reduces 
the Cauchy problem for the dKP equation to the Cauchy problem for the function 
F{^,y,t) introduced in (1.8), which is regular beyond tc- 


• in section 3 we study the singularity formation in the solution to the dKP equation 
as done in [35] , |38| . We then show that the local structure of the dKP solution near 
the point of gradient catastrophe, in a suitable system of coordinates, is equivalent 
to the unfolding of an A 2 singularity. We derive the self-similar structure of the 
lip-shaped domain where the solution of the dKP equation becomes multivalued. 


• In section 4 we give a heuristic derivation of the shock front position near the 
critical point of the solution of the dKP equation, and study the solution of the 
dKP equation when dissipation is added (called the dissipative dKP equation). 
Using multiple-scale analysis, we show that in the limit of small dissipation and 
near the critical point of the dKP solution, the solution of the dissipative dKP 
equation converges to a Pearcey integral. 


• In section 5 we compare our analysis with detailed numerical simulations. 
Solutions for initial data with and without symmetry with respect to y —y 
are studied. It is shown that our numerical approach allows to continue dKP 
solutions to a second gradient catastrophe, well after the first catastrophe has 
occurred. We find no indication for blow-up of the solution to the transformed 
dKP equation. 
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2. Solution by characteristic transformation 


We consider the Cauchy problem for the dKPI equation 

{ut ^yyi 

u{x, y,t = 0) = uo{x, y), a:, y G M, t G M"*". 


( 2 . 1 ) 


Since we are interested mainly in local properties of the solution, we will assume that 
uo{x,y) is in the Schwartz class, namely it is smooth and decreases rapidly at infinity. 


Equation (2.1) can also be written in the evolutionary form 


Ut + uu^ = Uyy, (2.2) 

where df^f{x) = f{x')dx'. This has the form of a nonlocal conservation law 


u 


rtt + Vf = 0, f = UyBy, 

with and By unit vectors in the x and y directions. As a result. 


(2.3) 


u{x,y,t)dxdy = / Uo{x,y)dxdy. 


(2.4) 


Similarly, 


{u% + 2u^/3 - {d^ ^UyY + {2ud^ ^Uy) = 0, 


and hence the norm is also a conserved quantity: 


(2.5) 


M{t) = / u^{x,y,t)dxdy = / u^{x,y)dxdy. 


( 2 . 6 ) 


Since the left hand side of (2.1) is a total derivative, solutions have to satisfy the 
constraint 

/ Uyy{x,y,t)dx — 0, t > 0. 

Jr 

If the initial data does not satisfy such constraint, a low decays at infinity occurs for 
t > 0 even for initial data in the Schwartz class. This is a manifestation of the infinite 
speed of propagation in the dKP equation. For this reason we choose initial data such 
that 

/ uo{x,y)dx = 0 (2.7) 

JR 

for all so that the dynamical constraint is satisfied also at t = 0. After these 


preliminaries we transform the dKP equation using (1.8), to find an equation for 

F{i,y,t). 
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Proposition 2.1 The equations (1.8) give a solution to the dKP equation with smooth 
initial data uo{x,y) in implicit form, if the function F{f^,y,t) satisfies the equation 




( 2 . 8 ) 


with initial data 


F{x,y,D) =uo{x,y). (2.9) 

Proof Differentiating the second equation in ( |1.8 ) with respect to x, t and y we find 

e. = ^ (2.10) 




where we have defined A = 1 + tF^. Thus the derivatives of u with respect to the 
variables are 

= + ( 2 . 11 ) 

and 

z? z? 

( 2 . 12 ) 


X ^ , 


A 


Now the Hopf equation becomes 


n -L 


(2.13) 


which confirms that F is time-independent in this case. Differentiating (2.12) a second 
time, we find 


Uyy — 


Fy_ 

A 


ZJL 

A 


A 


1 

A 


yy 


til 

A 


after some manipulations. But this means that if u{x,y,t) satisfies the dKP equation 
(2.1), F{f,y,t) satisfies (2.8) with initial condition (|2.9|). □ 


We rewrite the equation (2.8) in the evolutionary form 

F, = apF„ + t(dpF„F^ - F^) 

where is the inverse of a derivation. We observe from the above equation that the 
nonlinear terms are multiplied by the time t and this show that for small times the 
nonlinear effects are damped. This observation qualitatively explains the fact that the 
function develops a singularity after u{x^y^t) becomes singular. 

For the remainder of this paper we will focus on solutions to the transformed 
equation (2.8). We observe that (2.8) also conserves the integrals over F and 
which we will use to check our numerics. Namely for n integer one has 

[ u^dxdy= [ F^Ad^dy= [ F^d^dy +—^ [ {F^I^d^dy = [ F^d^dy. 

9r2 Jr2 Jr2 .V. I 1 


n + 1 
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Figure 1. A typical sequence of wave breaking, as described by showing the 

lip-shaped domain inside which the wave overturn s. T he singularity first appears at 
the origin, then spreads rapidly in the direction (3.3). The scale of the lip is 
in the x-direction, and PP in the ^-direction. Full lines are solutions of A = 0 at 
t = 0.01,0.1, and 0.4, while the dashed line is the shock front, which has to be inserted 
in accordance with the shock condition (4.5), to be discussed in Sectionbelow. 


In particular, conservation of the norm (2.6) gives the constraint 


V,t)di,dy= / ul(x, y, 0)dxdp. 


(2.14) 


The transformation (1.8) has been constructed so as to unfold the overturned 


profile onto the initial condition in the case of a ^-independent initial condition. It 
is thus intuitive that if the overturning is modulated in the ^-direction, it is unfolded 
onto a function which shows no overturning, and having a weak dependence 

on y and t only. 


3. Overturning of the profile 


For generic initial data the solution of the dKP equation encounters a gradient 
catastrophe at points where the transformation ( |1.8| ) is not invertible [HS] 

= '^FtF^{^,y,t) = 0, (3.1) 


and consequently the gradients and Uy go to infinity, cf. (2.12). This is illustrated 
in Fig. [^for generic initial data, based on the local description to be developed below. 
The singular time F where the gradient catastrophe occurs first is the smallest t such 
that holds. Since for t < 4 the quantity A(^, y, t) has a definite sign in the ^ and 
y plane, A{Fy,tc) must be a zero as well as an extremum: A = = 0. Thus 
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the two-dimensional gradient catastrophe is characterized by the equations: 

^ + = 0 
F^^ = 0 

F^y = 0 (3.2) 

u(x,y,t) = F{^,y,t) 

X = tF{^,y,t) +^. 

The first three equations of ( |3.2[ ) determine the coordinates and of the 

singularity in transformed variables, taken as the origin in Fig. The x and u 

coordinates are recovered by substitution into the last two equations. One finds 
that {tFydx + dy)u{x, y, t) = Fy < oo, hence there is no gradient catastrophe in the 
transversal direction characterized by the vector field 

tFyBx + Gy, ( 3 . 3 ) 

see Fig. For generic initial conditions, the second derivatives of A will be nonzero 
at the gradient catastrophe: 


Vc, tc) 7 ^ 0 , Vc, tc) 7 ^ 0 F^yy{^c, Vc, F) 7 ^ 0 . 


(3.4) 


The conditions ( |3.2| ),(3.4) correspond to a cusp singularity in the notation of P, and 
will be found to describe the generic singularity for the dKP solution. The condition 
that F remains smooth, and thus the right hand side Fyy of (2.8) is finite, results in 
the additional constraints 


ff + tciKY 


= 0, F‘, = 0, f,; + = 0, 


ty 


yy y 


(3.5) 


where with a super-script we indicate the derivatives evaluated at the critical point. 


We now give a local description of the two-dimensional wave front u{x, y, t), based 
on expanding F{F y, t) near the gradient catastrophe described by ( |3.2[ ). Our numerical 
simulations confirm that F(^, y, t) remains smooth in the (^, y) plane not only near the 
first singularity, but well beyond. The region where the wave is multivalued has the 
typical lip shape also seen in the caustic surface of light waves near the cusp catastrophe 
[H]. We will show them to be self-similar with width in the horizontal direction 
and P-F in the transversal direction where t = t — tc, as done in [35]. The same scalings 
have been observed in [13| in the context of the 2-dimensional kinematic wave equation. 


In order to illustrate the way in which (1.8) unfolds the singularity, it is instructive 
to consider a family of exact solutions to (2.1) obtained in [36| : 


u(x,y,t) 


- 2ut 

Vi [ 4t 


(3.6) 







Shock formation 


11 


where B is an arbitrary function of one variable. The validity of (3.6) can be checked 


explicitly by substitution. Clearly (3.6) can be re-parameterized in the form 

, 2 ' 




i = 2^/tB(C - |j) + C, 


(3.7) 


which shows that F, as defined by (1.10), is 




A singularity in the dKP solution occurs when the second equation in (3.7) is no 


longer invertible; the first time this occurs is the critical time tc > 0, determined by 

1 


y/F = min 
CeR 


2B. 


In order to write (3.6) in terms of our function we use the double re- 

parameterisation: 


u{x,y,t) = F{Fy,t), x = tF{Fy,t) +F (3.8) 

12 2 

= ? = ViB(C-|) + C (3.9) 


We observe that has a singularity when the second transformation in (3.9) is 

no longer invertible, namely at a critical time for the transformed equation (2.8) 


(f 


4tc- 


It is also straightforward to check that at F{Fy:t) satisfies the constraints (3.5). 
Indeed one calculates directly from (3.9) that 


B 2VtB' +1 ^ y^B' 1 ^ ^ yB' 

2tl ViB' + l 4tl ViB' + F ^ 2tl{ViB' + iy 


so at the critical time one obtains the relations 


Ff 


yl_ pc 

Aty y 2ty 


which satisfy the first of the constraints in (3.5); the remaining constraints (3.5) are 
checked analogously. 
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3.1. Local analysis 

In order to study the solution near the gradient catastrophe we expand the generalized 


characteristic equation (1.8) in a Taylor series near tc, Xc, Vc, Uc and ^c- Part of the 
analysis below is already contained in [33], [HEj- Introducing variables relative to the 
singularity as 

X :=x-Xc i:=t-tc, y ■■= y - Vc, C C - Cc 
we have argued that x ~ and y ~ Since Au ~ it follows from the first 


equation of (1.8) that ^ ~ P/^. Thus to be consistent, we include all terms up to 

0(p/2): 


X = i(f' + t,Fn + ty(F; + t,F‘,) + t, F‘y + + ^F^y’^ 

t 


+ + Fp )«+ y'X'Af + f ))■ 


(3.10) 


This suggests introducing the shifted variables (using 4 = 


c = Fn^+ 


pc 

pc 


y 




X - r(F“ + t,F‘) - mp + t,F‘,) - to (py + + Ip„p 


F-tc—— y Frc——yt 


cx-rpc ^2^ '2"'' 


3 






€ pc 


T=- 

k 


P 

f + p 


= I AkP 


pc 


pc 

^yy 


k = 


6 


SO that in the variable (3.10) takes the form 

- c’ + TC = X + »(f", y\(\tp + f)). 


6 


(3.11) 


(3.12) 


Using the estimates ^ ~ y ~ and x ~ identified previously, the scaling 

X ^ AX 
t — X^t 
y X^y 


C x^f, 


(3.13) 
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Figure 2. The universal s-curve described by (3.14); for T > 0 the profile turns over 
to form a multivalued region. 


in the limit A ^ 0 reduces (3.12) to the universal s-curve 

-e + TC = X 


(3.14) 


shown in Fig. It is easy to conhrm that the function C(^, T) defined by (3.14) solves 

Ct + CCx = 0 (3.15) 


with initial condition C{X,T = 0) = {—X)^, completing our task of reducing (2.1) 
locally to the Hopf equation. A gradient catastrophe is encountered for A = 0, T = 0, 
and C = 0. 


Using the identities (3.5), we can now calculate the solution to the dKP equation 
(2.1), valid near the singularity. To leading order in the limit A ^ 0, it is consistent to 
expand u{x, y, t) to linear order in y: 

u{x,y,t) -Uc = F{^,y,t) - ~ Ff^+Ffy = C,{X,T) + Py, (3.16) 

with 




Fk ■ 


(3.17) 


Thus putting u = u{x,y,t) — Uc, from (3.14) we find the local profile to be an 
s-curve, which has the universal similarity form: 


-{u- Hyf + T{u-i)g)= X. 


(3.18) 


This is the central result of our theoretical analysis; the formula (3.18) is the unfolding 
of an A 2 singularity. It is a complete description of the self-similar behavior of the 
dKP solution near its singularity for generic initial data. In the ^-independent case, 
(3.18) coincides with the usual result (1.6). We now derive the form of this multivalued 
valued region in the x,|/-plane, shown previously in Fig. 
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3.2. Multivalued region 


As seen in Fig. the function f = C,{X,T), described by the cubic equation (|3.14|), 

9A '-' 

becomes multivalued for T > 0. From —— = 0 it follows that 

dC 


T = X\ or C = ±yr73. 


(3.19) 


SO that for X in the interval 


3\/3 


2 j,3/2 ^ X < ^ 2^3/2^ 


3\/3 


the function f{X,T) is multivalued. 


Reversing the coordinate transformations (|3.11|), we can write the hrst equation 

(3.20) 


(3.19) in the form 


t = ^tca + 71/2) , 


where we have introduced the constants 

“ = 13 = ^, 


7 = 


fpc 

jyy 

pc 


(3.21) 


Alternatively, (3.20) could also have been derived from (3.1), and expanding F in a 
power series around the singular point. 


The a:-coordinate of the boundary of overturning can be found from (3.10), which 
using (3.20) can be simplihed to yield 

X = —a 


5 ?+ fsf)+f(f'-«q 7 ) + 

ty (q - (p^y + tp^f + 


2~ yy^ ' yyy^ / ■ (3.22) 

Equations ( 3.20| ) and (3.22) describe a curve in the {x,y) plane, parameterized by 
An example was shown previously in Fig. [^for several values of i — t — tc, showing its 
characteristic “lip” shape [3j. 

The overturned region starts from the singular point and then expands, as seen 
in Fig. [TJ To understand the scaling of this expansion, we introduce the independent 
variables 

X - t(F^ - tliF^f) - F{F^y + 


= 


Ti =yt 


Ilf 


^--3/2 


(3.23) 


Then the lip described by equations (3.20) and (3.22) is reduced to the time- 
independent similarity form: 


^F«(s2 + 2/5FiS + 7Fi2) = i 

Xi = —a F 5iYi F 


(3.24) 
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Figure 3. The symmetric lip according to (3.26), ending in a cusp 


with the additional constants 




pc _ 


5. = trK 


yyy' 


(3.25) 


This demonstrates that the gradient catastrophe in the dKP equation has a 
universal spatial signature, parameterized by the constants q;,/ 5,7, hi, and 82 , all of 
which can be computed in terms of the initial data and its derivatives at the point of 
gradient catastrophe {xc,tc,yc)- The scalings introduced in (3.23) imply that the lip 
expands as in the propagation direction, and as in the transversal direction, 
as announced previously. A characteristic feature is the cusp at the corner of the lip. 
This is seen most easily for initial data which is even in y, for which the description 
simplihes considerably. All odd derivatives in y vanish, and we obtain 


I + 7 ^ 1 ^) = 1 

\ a Q 


(3.26) 


shown in Fig. Analyzing the neighborhood of the point s = 0, one finds directly 
that 

= ±|(2Fi)^/2 (3.27) 

which is a generic 3/2 cusp [T5] . 
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4. Dissipative dKP equation and shock solutions 


The solution (3.18) constructed in the previous subsection is unphysical for t > 0, in 


that it does not assign a unique value of u to every point x, y in the plane. In principle, 
one can construct an infinity of single-valued solutions from it, by choosing different 
points at which to jump from one branch to the other. For conservation laws in many 
space dimensions, physically motivated constraints, known as generalized Rankine- 
Hugoniot jump conditions, have been introduced. As result, a weak solution of the 
equations (usually called the inviscid shock), is singled out uniquely [HSl EH] • 

Another way to select a unique solution after the singularity, is to consider a 


dissipative version of (2.1) with a viscous term added to it, which keeps the solution 


regular at all times. In the limit of vanishing viscosity e these regular solutions are 


expected to converge to (3.18), with a particular jump condition being selected. In 


this case the shock is called the viscous shock. In the field of hyperbolic equations the 
problem of showing that the inviscid shock is equal to the viscous shock has generated 
a huge literature. We only mention some important references in one dimension 
mm, and many space dimensions [TH]. Below we give an heuristic derivation of 
the equivalence of the inviscid and viscous shock for the dKP equation. 

We consider the dissipative form of the dKP equation 

{ut + uu^ - €{u^^ + cyyy))^ = Uyy, u{x, y, t = 0, c)) = uo{x, y), (4.1) 

with c > 0, which satisfies 

\ d f 

SW+ l u^{x,y,t)dxdy =-e{ul +cul) <t). (4.2) 

/ Ot J^2 

For given e-independent initial data, the solution u{x,y,t,e) of the dissipative equation 


(4.1) is expected to be approximated as e ^ 0 and t < tchy the solution u{x,y,t) of 


the dKP equation (2.1). 


4 . 1 . Shock position 


On the other hand, (2.4) is still satisfied at finite e, so (smooth) solutions of (4.1) still 
conserve u in the limit e ^ 0. From the condition that u be satisfied across a shock, we 
can use the generalized Rankine-Hugoniot jump conditions as in [HS], which determines 
the shock position (see also 0 )- Namely, if is the normal velocity of the shock, one 
obtains 

Vn{ui - U2) = (f • n)^ - (f • 11)2 , (4.3) 

where n is the normal to the shock front, and indices 1 and 2 denote values in front 
and in the back of the shock, respectively. Assuming that the shock position is given 
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by the curve Xs{y,t), and using the flux f from (2.3), this yields 


Xs{Ui -U 2 ) = — 


u“ — U 2 , dxs 


H-— / Uydx, 

dy 


(4.4) 


where xij^ are x-values approaching the shock from the front and from behind, 
respectively. 

Now the singular contribution to u across the shock can be written in the form 

u = U2-\-{ui- U2)9 {x - Xs{y,t)), 

where 0{x) is the Heaviside function and Ui ^2 become functions only of y and t on the 
shock front x = Xs{y,t). Hence 

Oxg 

Uy = {U 2 )y + {ui - U 2 )y 0 {x - Xs{y,t)) - («! - U 2 )-^S{x - Xs{y,t)), 


and from (4.4) the jump condition at the shock finally becomes 

2 


X., = 


Ui + U 2 


dxs 

dy 


(4.5) 


Note that the shock speed in the x-direction is not only an average between u-values 
in front and in the back of the shock as for the Hopf equation, but on account of the 
right hand side of (2.1) an additional term arises. 


Since we have mapped (2.1) locally to the Hopf equation (3.15), standard theory 
[Hj tells us that the shock should be at X = 0, according to (3.11) the equation for the 
front becomes 

Q yyy^ J ' 


x,(s, t) = t(F- + t,F:) + + t,F;,) + 4 (f^p + tF;,f 

Ij. -3 1+ -3 -j 

-tc , f: \^ y - tc ^ y — Ff-^yt. 

o pc \2^ o c pc y s pc y 

This equation indeed satisfies (|4.5|) to leading order, since 


(4.6) 


Xsiy, t) = + tcFf + ytcFf^ + 


(4.7) 


and 


dXs 

dy 


\ 2 

1 = [tc (q + F^y) + 0(t)] = -tcFf - tcFf^y + 0(t), 


(4.8) 

having used (3.5). On the other hand, C± = ±\/T at X = 0, and so according to (3.16) 

“1 + ^2 ipc , o- 

—-— = + By. 


Combining the last three equations one can see that the approximate shock front (4.6) 
satisfies (4.5) to leading order. In Fig. [^we have plotted (4.6) as the dashed line. 
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4-2. Shock structure 


Having found the shock position, we now investigate the inner structure of the shock, 


in case a small amount of viscosity is present. This is achieved by mapping (4.1) onto 


Burgers’ equation [48], which in addition to (3.15) contains a dissipative contribution. 


We are looking for a solution u{x,y,t-,e) of the dissipative dKP equation near the 
gradient catastrophe {xc,yc,tc) of the (inviscid) dKP equation. To this end we use the 
ansatz 

u{x,y,t;e) = Uc +h{X,T;e)+ y^, (4.9) 

with X and T defined in ( |3.11 ). Using the same scalings as before, and balancing 
Ut oc eUxx, we are led to the multiscale expansion 


h{X,T;e) = \'sH{X,r;e) + 0{X‘^), 

O 

X = XX, T = xH, e = xh, y = X^y, 


(4.10) 


and find the following theorem: 


Theorem 4.1 Let u{x,y,t;e) = Uc + h{X,T]e) + yS be a solution of the dissipative 
dKP equation with X and T defined in l[3.11\ ). Suppose that for \t — Ul small the 
limit 

H{X,T]e) = limA-U(AT, X^I^T]X^e) 


A-^O 


exists and the function H{X,T;e) satisfies the asymptotic conditions 

H(X,T,e) = =F + 0{\X\-i), \X\ -> oo 


(4.11) 


for each fixed T G M. Then the function H{X,T]e) satisfies the Burgers equation 

Hr + HH;^ = aH;^;^, a = ^{l + c{KFff) (4.12) 


with k defined in (3.11) 


Proof Inserting (4.9) into the dissipative dKP equation one obtains 


iHr+HHx-^{l + {pFf)^) + f, +y{Ff - 


= H- 


rx 


dt \dy J . 


o o SX^Hyt ( 1 o 2 

oy oy \oy J oy^ 


k \ \^yJ ^y ^y 


(4.13) 


A’ 
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Using (3.11), the constraints (3.5), and the substitution (4.10) one arrives at the relation 


(Hr + Mv - I (1 + = O(A^), 

which in the limit A ^ 0 shows that the derivative of (|4.12|) is equal to zero. In order 


to hx the integration constant we use the asymptotic condition (4.11). 


□ 


We remark that the asymptotic condition (4.11) implies that the local solution 


near the point of singularity formation, matches the outer solution given by (3.14) and 


(3.16). We conclude that near the gradient catastrophe, up to the constant term Uc 


as 


well as a term linear in y, in a suitable co-ordinate system the solution to the dissipative 
dKP equation reduces to the solution of the one-dimensional Burgers equation. We 
will argue below that the particular solution to the Burgers equations relevant near the 


critical point, and described by the asymptotic form (4.11), also satisfies the equation 


X = HT-H^ + QaHHx - 


ixx- 


(4.14) 


4 .. 2 .I. Burgers equation. To find the local solution near the shock, let us recall the 
solution to Burgers’ equation 

Vt-\-vv^ = uv^^, (4.15) 


with initial data uo(x), where is a positive constant. An exact solution is obtained 
via the Cole-Hopf transformation mm to give the formula: 

f°° G{r],x,t) 

v{x,t,n) = —2udx\og / e dr), (4.16) 


where 


G{r],x,t) = / vo{s)ds + 


{x — rjf 

2 r~ 


(4.17) 


For i/ —>■ 0 the leading contributions to the integral come from the neighborhood 
of the critical points of G, namely 

d^G{r),x,t)\rj=^ = uo(0 - = 0- (4.18) 

Let us assume first that there is only one such critical point, which means that using 
the method of the steepest descent [IH] , the integral can be approximated as 


/ 


e 2r. dr] 


Attu 


djGi^,x,t) 


e ^ , 


where ^ = ^{x,t) is a solution of (4.18). Direct evaluation of (4.16), using the 

+ Oiu^), 


characteristic condition (4.18), then yields the solution 


v{x,t,p) = Uo(0 + ^ 


+ 1 )^ 


(4.19) 
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whose leading order contribution in the limit —?• 0 is the solution of the Hopf equation 


by characteristics. In addition, (4.19), contains a linear correction coming from the 


viscosity. Alternatively, the term linear in n can also be obtained using perturbation 
theory. 


The approximation (4.19) remains valid as long as the function G{r),x,t) has 
an isolated generic critical point, before the appearance of a gradient catastrophe. 
However, after the critical triple point {xc, tf) of the Hopf equation, where no(^c)A + l = 
0 and Vq(^c) = 0, ( |4.18| ) has three solutions, as illustrated on the left of Fig. below. 
Near this point G{ri, x, t) can be expanded in a Taylor series as 




AG := G{ri,x,t) - G{^c,Cc,tc) ~ Vq{Q'-^ - ^ - vff) + vlt, 


4! 


2 tl 


tr 


where x = x — Xc, t = t — G, rj = rj — ^c, Vc = no(^c) and Vq(^c) > 0. Thus near such 
critical point the solution of Burgers’ equation can in the limit i/ ^ 0 be approximated 
by 


v(x, t, z/) ~ Uc — 2 pdx log / exp 


'"(r 

2 tl t 


^ - T\ 

X — Vet) 


dr). (4.20) 


Some rescaling leads to the following (see also [ 12 ]) 


Theorem 4.2 JH]/ Near a gradient catastrophe (xc, tc) for the solution of the Hopf 
equation Vt + vv^ = 0, the solution v{x,t, of (4-l^ admits the following expansion 


/ , /l/\V4 

n(x,t,i/) = Uc + (^-j U 


x — Vrt 


1 ^'' V 


+ 0(r^'/'). 


(4.21) 


where Vc = v{xcjtc), n = the funetion U=U{a^h) is defined by 

r*+oo 


U{a, h) = -2da log / 


(4.22) 


Remark 4.3 The function U{a,b) satisfies both the Burgers equation 

U, + UUa = Uaa 

and the non-linear ODE in the a-variable, containing 6 as a parameter [S] 

a = Uh-U^ + 6UUa - HJaa- (4.23) 

The behavior of U(a, b) is illustrated in Fig. for negative, positive, and vanishing 
values of reduced time b, performing the integral in (|4.22) numerically. For large 


o 
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Figure 4. The Pearcey function U (a, b), for three different values of b. 


and fixed b the integral (4.22) behaves as the root of the cubic equation ( 3.14| ) (see 
below) 

U{a,b) = T + kl 

O 


The integral in (4.22) is related to the standard Pearcey function m, which 
describes the diffraction pattern near a cusp caustic mi, by a complex rotation. The 
relation (4.23) is convenient in deducing the asymptotic properties of [/(a, 6); it follows 
from 

r*oo 7 

__e“ I “2^^6+420)^^ _ 

dz 


f 


In catastrophe theory [g the potential 

A(z) = - 2 z% + 4zo 


(4.24) 


(the weight in the exponent of (4.22)) is the standard unfolding of the cusp catastrophe, 
which is a co-dimension 2 singularity. For 6 < 0 (before the gradient catastrophe), there 
is only one critical point 

0 = ^ = 4^^ - 4^6 + 4a, (4.25) 

dz 

which is the case we considered before (see Fig. [^. Evaluating A at the critical point 

A = -3z^ + 2z%, a = -z^ + zb. (4.26) 


(4.25) yields 
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Figure 5. The critical contributions to the integral (4.22) near a cusp catastrophe, 


at constant reduced time b. On the left, the critical points; there is a unique solution 
for 6 < 0, and three solutions for |a| < 2(6/3)^/^ if 6 > 0. On the right, the argument 
A of the exponential; for 6 < 0 there is a single contribution, for 6 > 0 there are three 
contributions to a given value of a. 


For 6 < 0 this gives the single-valued curve shown on the right of Fig. which leads 


to the solution (4.19). 


If on the other hand 6 > 0 (after the gradient catastrophe), in the range 
|o| < 2(6/3)^/^ there are three critical points. Thus the integral (4.22) has three 


contributions, with different values of A (cf. ( 4.26| )), which he on a swallowtail figure, 
as shown on the right of Fig. The integral is dominated by the smallest value of A, as 
long as the solutions are well separated. This means we must have 6^1 (cf. Fig. |^, 
or ^ 1. Closer to gradient catastrophe, a more sophisticated asymptotics is 

needed, or one has to evaluate the integral numerically, as we will do below. However, 
outside of the region b < 1, the integral is dominated by either solution zi or Z 3 . The 
changeover occurs for a = 0 , where A(. 2 i) = A(. 23 ), namely on the line x — Vct = 0 . 
This is exactly the shock front near the gradient catastrophe (xc,tc). 


3 . 


4.1 


4 . 2 . 2 . Pearcey integral and dissipative dKP equation. Choosing A = e 4 in Theorem 
we obtain that the solution to the dissipative dKP equation satisfies in the rescaled 


variables (4.10) the Burgers equation (4.12) with e = 1. Furthermore for t < t^ such 


solution is asymptotic to the Hopf solution (3.14). Combining these observations with 


Theorem 4.2 and remark 4.3, we come up the following conjecture. 


Conjecture 1 Let us consider the double scaling limit e ^ 0, x ^ Xc, y ^ Vc o,nd 
t ^ tc in such a way that the ratios 
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remain bounded with X and T defined in (3.11). Then the solution u{x,y,t;e) of the 
dissipative dKP equation near the first singularity for the solution of the dKP equation 
is described by the expansion 


u{x, y, e) ~ Me + (j 


1/4J, (. 


X T 


where 


6 (l + c {tcFfY 


+ y/3 + 


a = e- 


pc J.4 


and the funetion U(a,b) is the Peareey integral defined in (4.22). 


(4.27) 


For y-symmetric initial data the expression (4.27) rednees to the form 

u{x, y, t; e) ~ Me + (7 / U [ - ^^ 3/4 -,-- 1 + 0(e / ), (4.28) 


with k defined in (3.11). The center of the (smooth) shock front is located at X = 0, 
as found previously in the inviscid limit. 


5. Numerical solution 


In this section we present numerical solutions of the transformed version (2.8) of 


the dKP equation, which remain smooth well beyond the gradient catastrophe of 


the original equation (2.1), as we will demonstrate below. In addition, we treat 


the dissipative dKP equation (4.1), whose solutions are also observed to remain 


smooth. We use a Fourier method for the spatial dependence, and an exponential 
time differeneing (ETD) scheme for the time dependence, as previously for the dKP 
equation [25] . 

Both equations are written in evolutionary form 


F, = apF„ + t(F^dpFy,-F‘), 


and 


Uf T UUx ^yy 4” ^ (^xx 4~ ^^yy) i 


(5.1) 


(5.2) 


with a small dissipation parameter e. In Fourier space, the antiderivatives and 
are represented as Fourier multipliers —i/k^ and —i/kx, respectively. Here /c^, k^, ky 
are the dual Fourier variables of x, y respectively, and the Fourier transform of a 
variable will be denoted by a hat. Thus (|5.1|) and (5.2) can be written in the form 


Ut = Fu + M{u), 


( 5 . 3 ) 
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where £ is a linear, diagonal operator, which is iky/k^ for (5.1), and iky/k^ — ekf. for 
(5.2), and J\f{u) is a nonlinear term. The idea of the ETD scheme to be used here is 
to treat the linear part of (5.3) exactly. We use the fourth order EDT method by Cox 
and Matthews [10] , but other schemes offer a very similar performance 


To satisfy the constraint (2.7) on the initial condition, we choose initial data as the 
derivative of a function from the Schwarz space of rapidly decreasing smooth functions. 
This is well suited to a Eourier method, since a Schwarz function can be continued as 
a smooth periodic function to within our finite numerical precision. However, the 
nonlocality of (5.1) and (5.2) implies that solutions will develop tails with an algebraic 
decrease towards infinity. This follows already from the Green function of the linearized 
equations (22] • It was shown in [211 El] that discontinuities at the boundaries of the 
computational domain can nevertheless be avoided by choosing a large enough domain, 
and one can achieve spectral accuracy (an exponential decrease of the numerical error 
with the number of Eourier modes) over the time scales considered. 


The antiderivative in both (5.1) and (5.2) leads to Eourier multipliers which are 
singular in the limit of small wave numbers. These terms are regularized in Fourier 
space by adding a term of the order of the machine precision ( 10 ^®here). In |26|. the 

dKP equation (2.1) was solved for which is possible since solutions maintain the 

property of being the derivative of a Schwarz function. Together with an exponential 
integrator treating the term iky/k^ explicitly, this addressed all numerical problems 
stemming from this singular operator. 


However, an explicit treatment of all singular terms is not possible for (5.1), since 
J\f is singular as well, which leads to numerical problems for k^ 0. This can be 
addressed by applying a Krasny filter [2H] : all Fourier coefficients with modulus smaller 


than some threshold (typically 10“^°) will be put equal to 0. In all cases considered, our 
numerical algorithm could now be continued well beyond the first gradient catastrophe. 
For longer times, the above mentioned algebraic tails will lead to a slower decrease of 
the Fourier coefficients and thus to numerical problems once the numerical errors are 
of the order of the Krasny filter. For long time computations, which are beyond the 
scope of the current paper, one would have to use considerably larger domains and 
higher resolutions, or alternatively a spectral approach as in [I]. 


The accuracy of the numerical solution to (2.8) was monitored via the decrease of 
the Fourier coefficients, and checking the conservation of the norm (cf. (2.6),(2.14)). 
To this end we compute 

M{t) 


S{t) = 1 - 


(5.4) 


M(0)’ 

whose time dependence will be a measure of the numerical error. As shown in 
the maximum error in F may well be one to two orders of magnitude greater than 
S, but within these limits S is nevertheless a reliable indicator of the accuracy, if the 
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Breaking event 

Initial data 

tc 

Xc 

Vc 

Uc 


First 

—6dx sech^A/a:^ + y^ 

0.222 

1.79 

0 

2.543 

1.227 

Second 

—Qdx sech^y/x'^ + 

0.300 

-2.033 

0 

-2.48 

-1.289 


Table 1. Critical parameters for 


initial data (5.5). 


Ulie iilSL LWU WclVe UieclKilig 


Fourier coefficients decrease sufficiently rapidly. 


5.1. Shock formation for symmetric initial data 

We begin with the simplest case of initial data symmetric with respect to y —>■ —y. We 
choose the same initial condition as [2H], 

(5.5) 


UQ{x,y) = 


who solved the dKP equation ( |2.1[ ) in its original form. Near the gradient catastrophe, 
( 2 . 1 ) develops a discontinuity, and the numerical scheme employed in [25] breaks down. 
By contrast, using the transformed equation (I^ 


we are able to reach the gradient 
catastrophe with much lower resolution (using serial instead of parallel computers), but 
are also able to continue the computation beyond the first and even secondary wave¬ 
breaking events. Beyond the gradient catastrophe, we identify the lines A = 0 along 
which the gradient of the solution blows up (cf. Fig. [^, and show that the solution 
of (2.8) yields the expected weak solution of dKP inside the lip region. We also show 


that the solution of ( 2 . 8 ) stays regular on time scales of order unity. 


In [25|, the first wave breaking event was observed at the critical time 4 = 
0.2216..., see Table Here we can identify tc directly from a solution of (2.8) by 
tracing the minimum of A over space. The first time this quantity vanishes or becomes 
just negative will be taken as the time tc. We use = Ny = 2 ® Fourier modes for 
x,y £ [— 57 r, 57 r]^ and Nt = 1000 time steps for t < 0.23. The first negative value of 
A is recorded for t = 0.222 ..., which is in agreement with [25] to within the accuracy 
of at least two digits. However, the present calculation requires much lower resolution 
to reach similar accuracy {N^ = Ny = 2^ compared to = Ny = 2 ^® in [25]), and 
accuracy can easily be improved. For example, after determining the critical time to 
a certain accuracy, one uses the required resolution in time close to the previously 
determined 4- This allows to determine the critical time with the same precision as 


the solution to (2.8), i.e., with the accuracy of the Krasny filter chosen here to be equal 
to 10“^°. For our purposes an accuracy of the order of 10“^ will be sufficient. 

The location of the critical point was identified in [25] as Xc = 1.79 ... and yc = 0. 
Here it is calculated for t = tc hj first finding the minimum ^c = 1.227. .., yc = 0 oi 
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t 



Figure 6. Measures of the smoothness of the solution to (2.8) with initial data (5.5). 
On the left, the time dependence of the maximum norm of F, as well as of and 
Fy; all decay for long times. On the right, the Fourier coefficients of the solution for 
t = 0.32. 


A, where Vc, 4) = 2.543 .... Then, using d^ , we find Xc = Cc + Vc, 4) = 

1.792..., again in excellent agreement with our previous result [2S], estimated to be 
correct to at least two digits. 

However, the solution F of ( |2.8| ) stays perfectly regular well beyond the critical 
time tc of the dKP solution u{x, y, t), as seen in Fig. On the left, we show that the 
maximum norms of the first derivatives of F remain bounded and smooth at tc, and even 
decay for long times (of course, the derivatives of the original variable u{x, y, t) diverge 
at a gradient catastrophe). On the right, for t = 0.32 we demonstrate exponential decay 
of the Fourier coefficients to the level of the Krasny filter, as expected for a smooth 
function. The relative norm 5{t) (cf. (5.4)) is conserved to the order of 10“^^. 
On account of the algebraic decay of the solution in Fourier space, the computation 
cannot be run for much longer than t = 0.35 at the current resolution. To be able to 
do so using a Fourier method, larger domains and higher resolution would be needed. 
However, there is no indication that the solution of (2.8) itself develops a singularity. 

Thus it is possible to continue the computation beyond the first wave breaking 
event, and to identify the second event, which occurs for negative x. This is of course 


not possible in the case of direct integration of (2.1) as in [25], where the numerical 
method fails at the first wave breaking. We use = 2^, Ny = 2^^ Fourier modes and 
W = 5000 time steps for t < 0.32. Proceeding as for the first break-up in tracing the 
minimum of A(^, y, t), we find F = 0.300 ... and x^ = —2.033 ..., see Table 

The corresponding profile u{x, y, t) can be seen in Fig. on the left. It is obtained 
by plotting F{Fy,t) (shown on the right) as a function of x = ^ -|- tF{Fy,t), as 
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Figure 7. Profiles obtained from a solution of the transformed equation (2.8) at 
t = 0.300, time of the second wave breaking event. On the left, the original solution 
for initial data (5.5); on the right, the profile u{x,y,t) obtained using the 
transformation ( |1.8[ ). The slices along the plane y = 0 (bottom) make it clear that 
the profile u{x,y,t) has overturned near x = 2 (first breaking), and is at the point 
of breaking near x = —2 (second breaking). The profile of F{^,y^t) remains smooth 
and single valued. 


required by ( |1.8D . For t > tc in di neighborhood of the blow-up point, one has that 
X = ^ + tF{^^y^t) is not invertible as a function of However we can still 

perform a parametric plot of u{x^y^t)^ which becomes a multivalued function in the 
region near the first critical point (xc,0,tc)- This is even clearer from the cut along 
the y = 0-axis shown on the bottom (recall that the critical points are all on the x- 
axis since the initial data are symmetric with respect to y ^ —y^ and since the dKP 
equation preserves this symmetry). Thus as for the solution to the Hopf equation via 
the characteristic method, a nonphysical solution which has overturned is obtained in 
the shock region. It is clear from the corresponding cut through F{^^y^t) shown on 
the bottom left that F remains smooth and single valued. 
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4.5 



’°f.82 1.825 1.83 1.835 1.84 1.845 1.85 
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Figure 8. On the left, the solution u{x, t) (blue line) of the dKP equation and 
its approximation (green line) ( |3.18| ) for t = 0.24 > tc ^ 0.222. The regions of 
multivaluedness of the solutions are projected on the (x, y)-plane. On the right 
top: a cut through u{x, y, t) along y = 0. On the right bottom: the corresponding 
multivalued regions of u{x, y) in the (x, j/)-plane (blue line: numerical solution; green 
line: local approximation.) 


We can now test to which extent the asymptotic description of the overturned 
region in Section]^ which only becomes exact in the limit t ~ tc; can approximate our 
numerical results. Recall that the profile is described by (3.18), while the shape of the 
overturned region is given by (3.23),(3.24). In Fig. we show a comparison between 
a numerical solution of the dKP equation, obtained through the transformation (1.8) 
(blue), with the local approximate solution (3.18) shown in green. At t = 0.24, i.e. 
shortly after overturning at tc = 0.222, there is good agreement in the description of 
the multivalued region. On the left, u(x, yt) is shown in a perspective plot, on the top 
right an s-curve is produced by a cut along the y = 0 plane. If corresponding cuts are 
considered for each value of r/, a lip-shaped region is obtained inside which the profile 
has overturned (bottom right). 

To test for the self-similar properties of the multivalued region, in Fig. [^we show 
the numerical result as function of the rescaled coordinates W, li, which are defined 
by ( |3.23| ) (red lines). Good agreement is seen with the asymptotic prediction (3.24) 
(blue fines), in particular for small values of three time distance f from the gradient 
catastrophe, as expected. The fact that the numerical results stay time independent to 
a good approximation demonstrates that the typical scales of the solution agree with 
the prediction (3.23): the width of the region scales like its height like . 
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Figure 9. Multivalued region of the solution of the dKP equation as found from 
= 0 for the initial data (5.5). Results are written in selfsimilar rescaled 
coordinates Xi and 15 defined by (3.23) for several values of i (red lines). The 
corresponding asymptotic boundary (3.24), shown in blue, is time-independent by 
construction. 



Figure 10. Numerical solution to the dissipative dKP equation (5.2) with c = 0 and 
e = 0.01 for initial data (5.5) at time t = 0.32 on the left, and the corresponding 
Fourier coefficients on the right. 


We now turn to the numerical solution of the dissipative dKP equation (4.1), and 
to the comparison with our asymptotic theory, which is given by (4.27) in the general 
case, and by (4.28) for symmetric initial data. To resolve the strong gradients in 
the solutions to the dissipative dKP equation (5.2) that occur for small e, much higher 
resolution is needed than for the solution of (5.1) for the same initial data. For e = 0.01 
(with c = 0) we use W = 2^^, Ny = 2^^ and W = 5000 to find the solution of (5.2) 
with initial data (5.5) at t = 0.32, shown in Fig. 10 on the left. At this value of e, the 
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Figure 11. Top left: numerical solutions to the dKP equation (blue) and to the 
dissipative dKP equation ( |5.2| ) (green), for c = 0 and e = 0.01, using symmetric 
initial data (5.5). Shown is a slice along the line y = 0 at t = 0.32 > tc = 0.222. 
Top right: the asymptotic approximations (3.16) and (4.27) to the same solutions; 
the dashed line marks the shock position X = 0. Bottom: the dotted lines mark 
the multivalued regions for t = 0.32, according to the numerical solution to the dKP 
equation (blue), and according to the asymptotic theory ( 3.26| ) (green). The green 
solid line is the asymptotic prediction for the shock front, as given by (4.6), and 
the blue solid line is a numerical estimate based on the inflection point of the dKP 
solution. 


total loss of the norm (cf. (4.2)) is of the order of 2%. A comparison between the 
dKP solution and the Fourier coefficients, shown on the right, decay to below 10“^°, 


as for the solutions to (5.1). To achieve higher resolutions, parallel computation would 
be needed. 


In Fig. (top left), we show a slice through the same dissipative solution at 
y = 0 (green line), together with the corresponding dKP solution, which has become 
multivalued, as i ~ 0.1. The dissipative solution exhibits a sharp front close to where 
the shock discontinuity is expected to be. Both curves are to be compared to our 


asymptotic results, shown on the top right, with the s-curve (3.18) shown in blue, and 


the dissipative asymptotics (4.28) in green. The sharp front is seen to be localized 
around the theoretical shock position, shown as the vertical dashed line. Since i 
is only moderately small, there exists a 30% difference in the height of the s-curve. 
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Figure 12. On the left, in blue the solution to the dissipative dKP equation (4.1) 
for e = 0.01 and the symmetric initial data (5.5) at the critical time tc = 0.222 
and near the critical point, and in green the asymptotic solution (4.27) given by the 
Pearcey integral. On the right the same plot along the line y = 0- The dashed blue 
line is the solution of dKP equation and the green dashed line is the solution of the 
approximation (3.16) to the dKP solution. 


but otherwise the overturning of the dKP equation is well reproduced. Within these 
limitations, the shape and width of the shock front, as well as the front position within 
the s-curve, are very well reproduced. 


In the bottom graph of Fig. we report the multivalued regions, as well as the 
position of the shock front, as given by the numerical solution (blue curves, with the 
shock front as the solid line), and our asymptotic theory (green curves, shock front 
solid). Once more, there is fair agreement in the shape and size of the lip-shaped 
multivalued regions (dashed lines), described by the dKP equation. The numerical 
shock position is estimated from the inflection point of the dKP solution, the theoretical 
prediction is the curve X = 0. 

In Fig. 


12 


we show the solution to the dissipative dKP equation (4.1) for e = 0.01 
and the asymptotic description (4.27) for the symmetric initial data (5.5) at the critical 
time in the vicinity of the critical point. While the asymptotic formula provides the 
best local approximation being best near the critical point, it can be seen to also 
correctly reproduce the ^/-dependence. 


The approximation is also valid for small, nonzero values of i as can be seen in 
Fig. where the same situation as in Fig. [^is shown on the slice y = 0 for several 
values of f. 
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t =-0.022 t =-0.000 
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Figure 13. Solution to the dissipative dKP equation (4.1) for e = 0.01 and the 
symmetric initial data (5.5) in blue, the Pearcey asymptotic solution (4.27) in green 
and the (weak) dKP solution dashed on the line y = 0 for several values of t. 


5.2. Nonsymmetric initial data 


In this section we consider two different initial profiles which are not symmetric with 
respect to y ^ —y. The first, 


tip,!/, 0) = ea* |(a:+1)(!/- l)e “’j, 


(5.6) 


still retains a radial symmetry for ^ oo. As seen in Table 5^, we can follow 

the evolution through two successive gradient catastrophes. The second profile, 


u 


(x, t/, 0) = 


(5.7) 


does not possess radial symmetry for large and we are able to compute the 


first gradient catastrophe only, whose critical parameters are also given in Table 5.2 


To solve the Cauchy problem with initial data (5.6) for the dKP equation (2.8), 


we use Nx = 2® and Ny = 2^^ Fourier modes for (x,|/) G [—57r,57r]^ and Nt = 5000 
time steps for t < 0.15. The first critical time is reached at 4 = 0.08323 ..., the second 


critical time is tc = 0.1070...; all other critical parameters are reported in Table 5^ 
The relative computed Lf norm is conserved to the order of 10“^^, and the Fourier 


coefficients decrease to the order of the Krasny filter as can be seen in Fig. 14 (left) 
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Breaking events 

Initial data 

tc 

Xc 

Vc 

Uc 


First breaking 

60 :^ 1 (x -b 1)(?/ - 1 

0.0832 

-1.210 

-0.368 

-4.958 

-0.798 

Second breaking 


0.1070 

2.004 

-0.368 

4.4066 

1.534 

First Breaking 


0.086 

0.088 

-0.245 

-1.477 

0.215 


Table 2. Critical parameters for the first two wave breaking events, with weakly 
asymmetric initial data (5.6). For the strongly asymmetric initial data (5.7) only the 
first breaking could be computed. 




Figure 14. Same as Fig.|^ but with initial data ( |5.6| ) (left). The Fourier coefficients 
on the right are shown for t = 0.15. 


As seen in the same figure on the left, the norm of the solution F and the norm 
of its gradient also appear to decrease for large t, so again there is no indication of a 
blow-up of the solution. However, to be able to run the code for longer times, larger 
computational domains would have to be used. 


On the left of Fig. 15, we trace the boundary of the multivalued regions of u{x, y, t) 
at four times shortly after the first gradient catastrophe; the times t relative to the 
singularity are reported on the top of each graph. On the right of the same figure, 
the same multivalued regions are plotted as functions of the rescaled coordinates Xi 
and Yi defined in (3.23). Once more, in the rescaled coordinates the shape of the 
multivalued region is almost constant, and agrees well with the theoretical prediction, 
shown in blue. Note the slight asymmetry of the lip shape with respect to the reflection 
symmetry y ^ —y. 

For the initial data (5.7), the code is run with = Ny = 2^^ Fourier modes on 
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Figure 15. Left: boundary of the multivalued region found from a numerical solution 
to the dKP equation for the initial data (5.6), for several values of t > tc = 0.08323 ... 
in the original (x,y) variables. Right: The red boundaries on the right are the same 
data represented in self-similar variables Xi and Yi as defined in (3.23), predicted 
to be time-independent by our asymptotic theory. The corresponding self-similar 
boundary, given by (3.24), is plotted in blue. 




Figure 16. Left: numerical solution to the dKP equation (2.1) for strongly 
asymmetric initial data (5.7) at t = 0.15 > tc = 0.087. Right: The corresponding 
contour of the multivalued region A(^,^,t) = 0 (red), compared to the asymptotic 
theory (3.24) (blue); the dashed line corresponds to X = 0 as given by (4.6). 


the same spatial domain as before, using Nt = 2000 time steps for t < 0.15. The first 
gradient catastrophe is found at tc = 0.087..., see Table for the remaining critical 
parameters. The solution at the final time (cf. Fig. left) is strongly asymmetric. 
This also implies an asymmetry of the tails of the solution and thus a stronger effect of 
the algebraic decay of the solution towards spatial infinity. The asymmetry of the tails 
of the solution also affects the Fourier coefficients. Despite a higher resolution than that 
of Fig. there are small contributions to the high wave number Fourier coefficients 
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Figure 17. Left: numerical solution to the dissipative dKP equation (5.2) with 
e = 0.04, c = 1, for initial data (5.7), at t = 0.15. Center: the corresponding Fourier 
coefficients. Right: a slice of the left plot along the line y = —0.4985 (green), together 
with the corresponding solution of the dKP equation (blue). 



Figure 18. In blue the solution of the dissipative dKP equation and in green the 
Pearcey asymptotic solution (4.27) for e = 0.01 and the strongly asymmetric initial 
data (5.7) at the critical time tc and near the critical point of the dKP solution. 


along the ky axis above the Krasny filter, which eventually cause the numerical scheme 
to break down. As a result, we do not reach a second catastrophe in this example. At 
t = 0.15, the relative computed LS norm is still conserved with an accuracy in the order 
of 10“^^. The norm of F and of its gradient do not indicate blow-up, but they are 
also not decreasing. If the solution exists for large t also, then the computation did 
not reach the asymptotic regime. 

The asymmetry of the solution can also clearly be seen from the contour delimiting 
the multivalued region, seen as the red line in Fig. (right). This is compared to the 
asymptotic theory at t = 0.063, shown as the blue line. Theory correctly describes the 
strong asymmetry and the orientation of the lip shape, but there are some quantitative 
differences. This indicates that the size of the critical region is smaller in the case of 
strong asymmetry. 
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For the dissipative dKP equation for the initial data (|5.7|), we consider e = 0.04 to 
obtain the solution shown in Fig. I 


I on the left. The Fourier coefficients in the middle 
of the same figure are also rather asymmetric, but decrease to the order of the Krasny 
filter. Due to the higher value of e, the loss of the norm is of the order of 22.2%. On 
the right of Fig. we compare the dissipative solution to the corresponding solution 
of the dKP equation. Although the width of the front is greater, owing to a higher 
value of e, it is set inside the s-curve where the shock position is expected to be. 


In Fig. 18 we show the dissipative dKP equation (4.1) for e = 0.01 for initial data 
(5.7). While in the symmetric case Fy = 0, here we have Fy ~ —17.39, consistent with 
a strongly asymmetric shock. Even in this case, the full two-dimensional structure of 
the step is well described by the asymptotic theory. 


6. Conclusions 


We have introduced a coordinate transformation, inspired by the method of 
characteristics, to investigate wave breaking in the dispersionless Kadomtsev- 
Petviashvili equation. As a result, the entire region where the profile is overturned is 
mapped onto a smooth and single valued function. The transformed equation remains 
smooth near the gradient catastrophe. Moreover, our numerics show that solutions 
remain smooth even beyond secondary wave breaking events. This permits us to 
compute solutions up to the first gradient catastrophe with much reduced numerical 
effort, and then to continue into the overturned region, where direct numerical 
simulations of the dKP equation fail. From the overturned profile, one can reconstruct 
the shock position, using the jump condition (4.5). 


Using the fact that the transformed profile remains smooth at the gradient 
catastrophe, we have calculated the local similarity form of the profile. This allows 
us to calculate the lip shape of the overturned region analytically, and to find the 
position of shock. Both the shape and the scaling properties of this region agree well 
with numerical simulations. 


We have also investigated the dissipative version of the dKP equation, which 
regularizes the gradient catastrophe. We performed direct numerical simulations of 
this equation for small dissipation, which we continued beyond the first gradient 
catastrophe. Results agree with expected shock solutions, except that the jump at 
the shock position is replaced by a smooth but rapidly varying profile. To investigate 
the shape of this profile, we use our characteristic transformation to map the dissipative 
KP equation locally to Burgers’ equation, which we can solve to obtain a local similarity 
description of the profile in two dimensions. Asymptotic analysis leads to a description 
of the profile in terms of Pearcey’s function, which is in good agreement with numerics. 
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We believe that the methods developed in this paper are of interest to study 
shock formation in a wider class of hyperbolic equations, including the compressible 
Euler equation. Here a significant complication lies in the fact that there are two 
families of characteristics in the corresponding one-dimensional problem, and hence a 
transformation based on a single characteristic cannot be expected to lead to a solution 
which avoids overturning for all times. However, shocks are generically expected to form 
with respect to one of the two characteristics only [30], so a transformation such as (1.8) 
will still be able to unfold the profile locally. However, the necessary transformation 
will depend on which of the characteristics is involved, and thus implicitly on initial 
conditions. 
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